Evaluation of Effectiveness of Global COVID-19 Vaccination Campaign

To model estimated deaths averted by COVID-19 vaccines, we used state-of-the-art mathematical modeling, likelihood-based inference, and reported COVID-19 death and vaccination data. We estimated that >1.5 million deaths were averted in 12 countries. Our model can help assess effectiveness of the vaccination program, which is crucial for curbing the COVID-19 pandemic.

where S, E, I, H, D, and R denote the mean number of the population that are susceptible (S), exposed (E), infectious (I), a delayed class of hospitalized (H) persons between infectious (I) and dead (D) and recovered (R), respectively. η is the proportion that become fully protected after vaccination (a proxy measure of vaccine efficacy). β is the transmission rate, σ is the rate at which exposed persons become infectious, γ is the recovery rate, π is the proportion of infectious persons that enter the delayed class H, κ the rate persons discharged from the delayed H class, and θ the proportion of deaths of the discharged persons. We suppose that reinfection and breakthrough infections are not significant in terms of contributing to deaths, which is a reasonable first approximation.
Following the literature, σ = 0.5 per day, γ = 0.33 per day and κ = 0.0833 per day, such that the mean generation time (sum of mean latent period and mean infectious period) is 5 days (4,5), and the mean duration from infection to death is 17 days, which are largely in line with observations (6). It is assumed that π = 0.15 for countries other than India where π = 0.03 (7,8), and and are estimated through fitting. The choice of π is not our immediate interest because we do not fit data for ; hence it is not dealt with in our fitting exercise. The infection fatality rate (IFR) = θπ. The choice of π = 0.03 in India was based on considerations of the infection to case ratio there (as high as 24:1), low reported deaths and high estimates of seroprevalence in India (8). Namely the reported deaths are relatively low per capita and all serologic studies in India suggest a large proportion of the population have been infected. We assume θ is in the range of 0.02-0.04, which means an IFR in the range of 0.3%-0.6% in the 11 countries except for India, which is reasonable (7). The IFR in India is 1/5 of that in the other 11 countries (8).
We allowed the time-dependent transmission rate, β(t), to be estimated by an exponential cubic spline (9) with several nodes nβ = 10 and an upper limit of 608, such that the reproduction number (without vaccine) is in a reasonable range with an upper limit of 5. The choice of cubic spline was the same as in our previous studies in modeling multiple waves of infections (10)(11)(12).
Alternatively, one could use the mobility index data in the transmission instead of cubic spline.
However, the mobility index data alone are insufficient (7,13,14). The emergence of new variants with increased transmissibility will increase the overall transmissibility in our 1-strain model.
Because the risk for infection is not uniform in the population and some persons might have strong protection compared with others, we assumed for initial conditions that 5% of the population were somehow protected or possibly had pre-existing cross-immunity from other coronaviruses (15). The initial E and I populations were equal and randomly chosen in the range of 0-10,000. The H class was given 1/10 the population of the I class, and the D class had 1/10 the population in the H class.
A partially observed Markov process (POMP) model with a maximum likelihood based iterated filtering technique was used to fit the mortality data (11). As mentioned, the transmission rate, β(t), was taken as an exponential cubic spline (9) to account for the simultaneous impact of all possible interventions excluding vaccination. The fitting procedure can be found at https://kingaa.github.io/sbied. Appendix Figure 1 shows the fitting and simulation results of model 1 with vaccine efficacy (VE) set at 85%.

Model 2
In model 2, we extended model 1 by including an additional vaccinated compartment (V) for tracking the dynamics of vaccinated but only partially susceptible persons (16). Thus, we further consider reduced susceptibility, reduced fatality rate due to vaccination, or both. The equations for model 2 are: Here, ψ is the parameter that accounts for the reduced susceptibility of vaccinated persons, where 0<ψ<1. We show fitting results of model 2 with ψ = 0.6 (Appendix Figure 2).

Vaccination Rate
We downloaded data for the vaccination rate, v(t), from the Our World in Data Web site (17,18), which is the proportion of the whole population vaccinated per unit of time. First, we calculated ṽ(t), the proportion of susceptible persons vaccinated per unit of time. The population is divided into 2 groups, vaccinated and unvaccinated; vaccination is only delivered to the unvaccinated group, which includes both susceptible and recovered persons. The rate at which susceptible persons are vaccinated is given as where is in units of days (1)(2)(3)16). We assume a delay of 14 days between the delivery of the second vaccine dose and the onset of protection, thus:

Asymptomatic Cases
A large proportion of infections are asymptomatic and less infectious than symptomatic cases, as reported in our earlier works (19,20). However, we adopted a simple homogeneous model that aggregates both the symptomatic and asymptomatic cases following other previous studies, such as Yang and Shaman (13).

Discussion
It is possible that population behavioral patterns become more careless and unstable due the widespread availability of vaccines over time (21). This might modulate the transmissibility across the epidemic and consequently cause us to overestimate the total deaths averted because of a vaccination campaign. To assess the effects of the vaccination, we compared the scenarios of with-vaccination (baseline scenario) and without-vaccination (counterfactual scenario). To test the sensitivity of varied transmissibility, we considered 5 sets of simulations all without-vaccination (v(t) = 0), but the transmission rate after April 16, 2021 was reduced by 0%, 10%, 15%, 20%, and 50% of the baseline scenario's level (reconstructed transmissibility from data with vaccination). We plotted the number of deaths that would have been averted as a percentage of the total population of each country with model 1 (VE = 85%) and these 5 counterfactual scenarios on transmissibility (Figure in main text).
Thus, if the reduction in model transmissibility is very large, say 50%, the disease will go extinct and few persons will die, which is not so different to the scenario under a successful vaccination policy. As such, we saw virtually no difference between the model simulations with We examined what happens in between those 2 extremes. Of note, according to Figure in main text, a 20% transmission reduction is not enough to bring about disease extinction and the I class is still able to grow exponentially in some phases for many countries. As such, we saw that for many countries, such as the United Kingdom, Spain, Germany, the United States, Italy, and France, this 20% reduction in transmission is not very different in terms of deaths averted than the 0% transmission reduction, and implies our overestimation not too large. We saw <15% difference in the deaths averted in 8 of the 12 countries, namely the United Kingdom, Italy, Russia, France, the United States, Spain, Germany, and Canada. On the other hand, with the 20% transmission reduction in 2 countries, Mexico and Columbia, the herd immunity threshold was crossed and the disease rapidly became extinct. This indicates that a 20% reduction in the transmission rate is probably too large to be reasonable, and that level of reckless behavior is unrealistic, which was confirmed by examining a scenario of 25% transmission reduction, which led to disease extinction in most countries.
Other than the above, we know of no other method to explore the effects of reckless behavior that might lead to overestimations but recognize this as a possible limitation of the method.

Sensitivity Analysis
In the above, we consider model 1 with VE = 85% and model 2 with susceptibility reduction ψ = 0.6. Here, we consider variations in the model. We consider model 1 with VE = 75% and VE = 95%. We consider model 2 with ψ = 0.8. In model 2, we further replace θ with the following equation: Namely, we assume that the death rate, θ, drops while the proportion of vaccinated persons increases at a rate of the following:  Table 1).
We fit these 6 model variations (including our baseline model, which is model 1 version 1 with VE = 85%) to the respective data to find the maximum-likelihood parameter estimates.
All model variations fit the data reasonably well (Appendix Figure 2). We compared the model- From this and a closer examination (Appendix Tables 2, 3), we concluded that our estimates of deaths averted show reasonable robustness to changes in the model structure and parameters. We have further confirmed this with a study of much larger number of model variations than reported here. *Model 1 is a Susceptible-Exposed-Infectious-Hospitalized-Died-Recovered (SEIHDR) model; model 2 is an extension of model 1 in which the vaccinated group (V) has reduced susceptibility (controlled by ψ) and reduced death rates (controlled by ε). η is the proportion of the population that becomes fully protected after vaccination, a proxy measure of vaccine efficacy.